1 Analysing the effect of treatment on molecular scores using aligned-rank transformed analysis of variance

1.1 INTRODUCTION

This analysis assesses how felzartamab therapy affects molecular sores measured with gene expression data from biopsies from kidney allografts taken from patients with antibody mediated rejection. The assumptions of normality and heterogeneity of variance differ across molecular scores, thus we avail to alinged-rank transformation to carry out a non-parametric ANOVA. This is powerful in this context because it allows us to test for the interactive effect of treatment (i.e., placebo vs felzartamab) and time (i.e., change in scores from biopsy to biopsy). Aligned-rank transformation also allows for us to specify a mixed-model to handle repeated measurements (i.e., biopsies) across patients.

References:

• Wobbrock J, F.L., Gergle D, Higgins J The Aligned Rank Transform for Nonparametric Factorial Analyses Using Only ANOVA Procedures. In Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ’11) 143–146. doi:10.1145/1978942.1978963, https://depts.washington.edu/acelab/proj/art/.(2011).

• Kay M, E.L., Higgins J, Wobbrock J. ARTool: Aligned Rank Transform for Nonparametric Factorial ANOVAs. doi:10.5281/zenodo.594511, R package version 0.11.1, https://github.com/mjskay/ARTool. (2021).

• Elkin, L.A., Kay, M., Higgins, J.J. & Wobbrock, J.O. An Aligned Rank Transform Procedure for Multifactor Contrast Tests. in The 34th Annual ACM Symposium on User Interface Software and Technology 754–768 (Association for Computing Machinery, Virtual Event, USA, 2021).

1.2 DATA WRANGLING


First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# Load necessary libraries
library(tidyverse)
library(broom)
library(rstatix)
library(flextable)
library(officer)
library(emmeans)
library(multcomp)
library(ARTool)
library(rcompanion)
library(Biobase)
library(ggpubr) 
library(patchwork) 
library(ggprism) 

# Custom operators
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0

# Suppress dplyr reframe info
options(dplyr.reframe.inform = FALSE)

# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/results/flextable_artANOVA.RData")
# load plot data
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")

# Seed for reproducibility
set.seed(42)


Now we define the categories of molecular scores for easier organization downstream.

# Feature categories
vars_cfDNA <- c("cfDNA_cpml")
vars_abmr <- c("ABMRpm", "ggt0", "ptcgt0", "DSAST", "NKB")
vars_tcmr <- c("TCMRt", "tgt1", "igt1", "QCAT", "TCB")
vars_injury <- c("IRRAT30", "IRITD3", "IRITD5", "cigt1", "ctgt1")
vars_parenchyma <- c("KT1", "KT2")
vars_macrophage <- c("AMAT1", "QCMAT")

# DEFINE VARIABLES TO ASSESS ####
vars <- c(vars_cfDNA, vars_abmr, vars_tcmr, vars_macrophage, vars_injury, vars_parenchyma)


Now we process the data so we can iteratively carry out the analyses on each score.

# DEFINE THE DATA ####
data <- set %>%
    pData()  %>% 
    dplyr::select(Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, all_of(vars)) %>%
    dplyr::filter(Patient %nin% c(15, 18)) %>%
    left_join(., summarise(., sample_pairs = n(), .by = Felzartamab_Followup), by = "Felzartamab_Followup") %>%
    relocate(sample_pairs, .after = Felzartamab_Followup) %>%
    arrange(Felzartamab, Patient, Group)


# WRANGLE THE PHENOTYPE DATA ####
data_00 <- data %>%
    expand_grid(category = c("cfDNA", "ABMR", "TCMR", "macrophage", "injury", "parenchyma")) %>%
    nest(.by = category) %>%
    mutate(
        features = map(
            category,
            function(category) {
                if (category == "cfDNA") {
                    vars_cfDNA
                } else if (category == "ABMR") {
                    vars_abmr
                } else if (category == "TCMR") {
                    vars_tcmr
                } else if (category == "macrophage") {
                    vars_macrophage
                } else if (category == "injury") {
                    vars_injury
                } else if (category == "parenchyma") {
                    vars_parenchyma
                } 
            }
        ),
        data = pmap(
            list(features, data),
            function(features, data) {
                data %>%
                    dplyr::select(
                        Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, sample_pairs,
                        all_of(features)
                    )
            }
        )
    )


# WRANGLE THE DATA FOR UNIVARIATE TESTS ####
data_01 <- data_00 %>%
    mutate(
        data_univariate = pmap(
            list(data, features),
            function(data, features) {
                data %>%
                    pivot_longer(
                        cols = all_of(features),
                        names_to = "variable",
                        values_to = "value"
                    ) %>%
                    nest(.by = variable) %>%
                    mutate(
                        variable = variable %>%
                            factor(
                                levels = vars
                            ),
                        annotation = case_when(
                            variable %in% vars_cfDNA ~ "cfDNA",
                            variable %in% vars_tcmr ~ "TCMR-related",
                            variable %in% vars_abmr ~ "ABMR-related",
                            variable %in% vars_macrophage ~ "macrophage-related",
                            variable %in% vars_injury ~ "injury-related",
                            variable %in% vars_parenchyma ~ "parenchyma-related",
                            TRUE ~ " "
                        ) %>%
                            factor(
                                levels = c(
                                    "cfDNA",
                                    "ABMR-related",
                                    "TCMR-related",
                                    "macrophage-related",
                                    "injury-related",
                                    "parenchyma-related",
                                    "archetypes",
                                    "rejection PC",
                                    "injury PC"
                                )
                            ),
                        score = case_when(
                            variable == "cfDNA_cpml" ~ "Donor-derived cell-free DNA (dd-cfDNA, cp/mL)",
                            variable == "TCMRt" ~ "TCMR classifier (TCMRProb)",
                            variable == "TCB" ~ "T cell burden (TCB)",
                            variable == "tgt1" ~ "Tubulitis classifier (t>1Prob)",
                            variable == "igt1" ~ "Interstitial infiltrate classifier (i>1Prob)",
                            variable == "QCAT" ~ "Cytotoxic T cell-associated (QCAT)",
                            variable == "ABMRpm" ~ "ABMR classifier (ABMRProb)",
                            variable == "DSAST" ~ "DSA-selective transcripts (DSAST)",
                            variable == "NKB" ~ "NK cell burden (NKB)",
                            variable == "ggt0" ~ "Glomerulitis classifier (g>0Prob)",
                            variable == "ptcgt0" ~ "Peritubular capillaritis classifier (ptc>0Prob)",
                            variable == "AMAT1" ~ "Alternatively activated macrophage (AMAT1)",
                            variable == "QCMAT" ~ "Constitutive macrophage (QCMAT)",
                            variable == "IRITD3" ~ "Injury-repair induced, day 3 (IRITD3)",
                            variable == "IRITD5" ~ "Injury-repair induced, day 5 (IRITD5)",
                            variable == "IRRAT30" ~ "Injury-repair associated (IRRAT30)",
                            variable == "cigt1" ~ "Fibrosis classifier (ci>1Prob)",
                            variable == "ctgt1" ~ "Atrophy classifier (ct>1Prob)",
                            variable == "KT1" ~ "Kidney parenchymal (KT1)",
                            variable == "KT2" ~ "Kidney parenchymal - no solute carriers (KT2)"
                        ), .before = 1
                    ) %>%
                    arrange(annotation, variable)
            }
        )
    ) %>%
    dplyr::select(category, data_univariate) %>%
    unnest(data_univariate) %>%
    mutate(n_cat = score %>% unique() %>% length(), .by = category, .after = variable)


This is what the data look like after we have organized them. Each row of the dataframe contains the data for each score, the category it belongs to, its category annotation, its long-format name, and the name of the variable in the data.

data_01 %>% print(n="all")
## # A tibble: 20 x 6
##    category   annotation         score         variable n_cat data    
##    <chr>      <fct>              <chr>         <fct>    <int> <list>  
##  1 cfDNA      cfDNA              Donor-derive~ cfDNA_c~     1 <tibble>
##  2 ABMR       ABMR-related       ABMR classif~ ABMRpm       5 <tibble>
##  3 ABMR       ABMR-related       Glomerulitis~ ggt0         5 <tibble>
##  4 ABMR       ABMR-related       Peritubular ~ ptcgt0       5 <tibble>
##  5 ABMR       ABMR-related       DSA-selectiv~ DSAST        5 <tibble>
##  6 ABMR       ABMR-related       NK cell burd~ NKB          5 <tibble>
##  7 TCMR       TCMR-related       TCMR classif~ TCMRt        5 <tibble>
##  8 TCMR       TCMR-related       Tubulitis cl~ tgt1         5 <tibble>
##  9 TCMR       TCMR-related       Interstitial~ igt1         5 <tibble>
## 10 TCMR       TCMR-related       Cytotoxic T ~ QCAT         5 <tibble>
## 11 TCMR       TCMR-related       T cell burde~ TCB          5 <tibble>
## 12 macrophage macrophage-related Alternativel~ AMAT1        2 <tibble>
## 13 macrophage macrophage-related Constitutive~ QCMAT        2 <tibble>
## 14 injury     injury-related     Injury-repai~ IRRAT30      5 <tibble>
## 15 injury     injury-related     Injury-repai~ IRITD3       5 <tibble>
## 16 injury     injury-related     Injury-repai~ IRITD5       5 <tibble>
## 17 injury     injury-related     Fibrosis cla~ cigt1        5 <tibble>
## 18 injury     injury-related     Atrophy clas~ ctgt1        5 <tibble>
## 19 parenchyma parenchyma-related Kidney paren~ KT1          2 <tibble>
## 20 parenchyma parenchyma-related Kidney paren~ KT2          2 <tibble>

1.3 CALCULATE MEDIAN SCORES AND MEDIAN EFFECT OF TREATMENT


At this stage we can summarize the median values for each score by treatment and follow-up.


(Note that the median scores are used for general interpretation of treatment effects because aligned-ranks, the data being assessed for inference regarding the effects of treatment, are not easily interpretable.)

# UNIVARIATE MEDIANS ####
data_02 <- data_01 %>%
    mutate(
        medians = map(
            data,
            function(data) {
                data %>%
                    reframe(
                        median = value %>% median(),
                        IQR = value %>% IQR(),
                        sample_pairs = n(),
                        .by = c(Followup, Felzartamab, Felzartamab_Followup)
                    ) %>%
                    arrange(Felzartamab_Followup)
            }
        ),
        medians_delta = map(
            medians,
            function(medians) {
                medians %>%
                    reframe(
                        median_delta = combn(median, 2, diff) %>% as.numeric(),
                        IQR_delta = combn(IQR, 2, mean) %>% as.numeric(),
                        sample_pairs = sample_pairs,
                        .by = c(Felzartamab)
                    ) %>%
                    mutate(
                        Followup_pairwise = rep(c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52"), 2) %>%
                            factor(levels = c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52")),
                        .before = sample_pairs
                    ) %>%
                    mutate(
                        median_delta_delta = combn(median_delta, 2, diff) %>% as.numeric(),
                        IQR_delta_delta = combn(IQR_delta, 2, mean) %>% as.numeric(),
                        .by = c(Followup_pairwise),
                        .before = sample_pairs
                    )
            }
        )
    )


Let’s take a look at what the median summaries look like for the ABMRprob score.

data_02 %>%
    dplyr::filter(variable == "ABMRpm") %>%
    pull(medians) %>%
    pluck(1) %>%
    flextable::flextable()

Followup

Felzartamab

Felzartamab_Followup

median

IQR

sample_pairs

Baseline

Placebo

Baseline_Placebo

0.6695

0.35525

10

Week24

Placebo

Week24_Placebo

0.7740

0.39025

10

Week52

Placebo

Week52_Placebo

0.5185

0.49825

10

Baseline

Felzartamab

Baseline_Felzartamab

0.7360

0.30775

10

Week24

Felzartamab

Week24_Felzartamab

0.1675

0.35700

10

Week52

Felzartamab

Week52_Felzartamab

0.7025

0.53875

10


Now we can look at the change in median ABMRprob scores in the placebo and felzartamab arms, as well as the differnce in those medians (i.e., ΔΔ median - this is our summarized treatment effect).

data_02 %>%
    dplyr::filter(variable == "ABMRpm") %>%
    pull(medians_delta) %>%
    pluck(1) %>%
    flextable::flextable()

Felzartamab

median_delta

IQR_delta

Followup_pairwise

median_delta_delta

IQR_delta_delta

sample_pairs

Placebo

0.1045

0.372750

Baseline - Week24

-0.6730

0.3525625

10

Placebo

-0.1510

0.426750

Baseline - Week52

0.1175

0.4250000

10

Placebo

-0.2555

0.444250

Week24 - Week52

0.7905

0.4460625

10

Felzartamab

-0.5685

0.332375

Baseline - Week24

-0.6730

0.3525625

10

Felzartamab

-0.0335

0.423250

Baseline - Week52

0.1175

0.4250000

10

Felzartamab

0.5350

0.447875

Week24 - Week52

0.7905

0.4460625

10

1.4 CARRY OUT ALIGNED-RANK TRANSFORM ANOVA


Now we carry out the analysis. This includes testing the specific contrasts for the interaction of treatment and time (i.e., our effect of treatment).

# UNIVARIATE NONPARAMETRIC TESTS ####
artANOVA <- data_02 %>%
    mutate(
        art = map(data, art, formula = value ~ Followup * Felzartamab + (1 | Patient)),
        art_aov = map(art, anova, type = "II"),
        art_aov_tidy = map(art_aov, tidy) %>% suppressWarnings(),
        art_aov_contrast = map(
            art,
            art.con,
            formula = "Followup:Felzartamab", adjust = "fdr", method = "pairwise", interaction = TRUE, response = "art"
        ),
        art_aov_contrast_tidy = map(art_aov_contrast, tidy),
        art_aov_contrast_cld = map(
            art_aov_contrast_tidy,
            function(art_aov_contrast_tidy) {
                art_aov_contrast_tidy %>%
                    as.data.frame() %>%
                    cldList(adj.p.value ~ Followup:Felzartamab, data = .)
            }
        )
    )

1.5 CREATE A SUMMARY TABLE OF ANOVA RESULTS FOR EACH SCORE

# CREATE FLEXTABLE SUMMARY OF ART MODELS ####
title_art <- paste("Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients")
artANOVA %>%
    dplyr::select(annotation, n_cat, score, art_aov) %>%
    unnest(everything()) %>%
    dplyr::rename(p.value = `Pr(>F)`) %>%
    mutate(FDR = p.value %>% p.adjust(method = "fdr"), .by = annotation) %>%
    dplyr::select(annotation, score, Term, `F`, p.value, FDR) %>%
    flextable::flextable() %>%
    flextable::add_header_row(values = rep(title_art, ncol_keys(.))) %>%
    flextable::merge_h(part = "header") %>%
    flextable::merge_v(j = 1:2) %>%
    flextable::fontsize(size = 8, part = "all") %>%
    flextable::align(align = "center", part = "all") %>%
    flextable::bg(bg = "white", part = "all") %>%
    flextable::bg(i = ~ FDR < 0.05 & Term == "Followup:Felzartamab", j = 3:6, bg = "#fbff00") %>%
    flextable::colformat_double(j = 2:3, digits = 2) %>%
    flextable::colformat_double(j = 4:ncol_keys(.), digits = 3) %>%
    flextable::border_remove() %>%
    flextable::bold(part = "header") %>%
    flextable::padding(padding = 0, part = "all") %>%
    flextable::border(border = officer::fp_border(), part = "all")  %>%
    flextable::autofit()

Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients

annotation

score

Term

F

p.value

FDR

cfDNA

Donor-derived cell-free DNA (dd-cfDNA, cp/mL)

Followup

4.418

0.019

0.029

Felzartamab

1.598

0.222

0.222

Followup:Felzartamab

10.396

0.000

0.001

ABMR-related

ABMR classifier (ABMRProb)

Followup

6.659

0.003

0.007

Felzartamab

1.620

0.219

0.235

Followup:Felzartamab

9.605

0.000

0.001

Glomerulitis classifier (g>0Prob)

Followup

11.768

0.000

0.001

Felzartamab

1.511

0.235

0.235

Followup:Felzartamab

10.743

0.000

0.001

Peritubular capillaritis classifier (ptc>0Prob)

Followup

14.188

0.000

0.000

Felzartamab

3.178

0.092

0.125

Followup:Felzartamab

7.749

0.002

0.004

DSA-selective transcripts (DSAST)

Followup

10.508

0.000

0.001

Felzartamab

1.747

0.203

0.234

Followup:Felzartamab

4.179

0.023

0.035

NK cell burden (NKB)

Followup

6.512

0.004

0.007

Felzartamab

1.956

0.179

0.224

Followup:Felzartamab

5.114

0.011

0.018

TCMR-related

TCMR classifier (TCMRProb)

Followup

3.555

0.039

0.165

Felzartamab

1.417

0.249

0.468

Followup:Felzartamab

1.124

0.336

0.504

Tubulitis classifier (t>1Prob)

Followup

1.016

0.372

0.506

Felzartamab

1.853

0.190

0.408

Followup:Felzartamab

0.065

0.937

0.968

Interstitial infiltrate classifier (i>1Prob)

Followup

3.414

0.044

0.165

Felzartamab

0.997

0.331

0.504

Followup:Felzartamab

0.032

0.968

0.968

Cytotoxic T cell-associated (QCAT)

Followup

5.693

0.007

0.107

Felzartamab

4.118

0.057

0.172

Followup:Felzartamab

0.927

0.405

0.506

T cell burden (TCB)

Followup

4.360

0.020

0.151

Felzartamab

2.868

0.108

0.269

Followup:Felzartamab

0.301

0.742

0.856

macrophage-related

Alternatively activated macrophage (AMAT1)

Followup

1.549

0.226

0.304

Felzartamab

1.047

0.320

0.320

Followup:Felzartamab

2.406

0.105

0.209

Constitutive macrophage (QCMAT)

Followup

3.378

0.045

0.209

Felzartamab

3.100

0.095

0.209

Followup:Felzartamab

1.426

0.253

0.304

injury-related

Injury-repair associated (IRRAT30)

Followup

0.055

0.946

0.946

Felzartamab

3.842

0.066

0.223

Followup:Felzartamab

1.915

0.162

0.347

Injury-repair induced, day 3 (IRITD3)

Followup

0.174

0.841

0.946

Felzartamab

3.268

0.087

0.223

Followup:Felzartamab

3.125

0.056

0.223

Injury-repair induced, day 5 (IRITD5)

Followup

0.073

0.930

0.946

Felzartamab

1.088

0.311

0.583

Followup:Felzartamab

2.586

0.089

0.223

Fibrosis classifier (ci>1Prob)

Followup

0.360

0.700

0.876

Felzartamab

3.482

0.078

0.223

Followup:Felzartamab

0.938

0.401

0.668

Atrophy classifier (ct>1Prob)

Followup

0.568

0.572

0.780

Felzartamab

3.297

0.086

0.223

Followup:Felzartamab

0.602

0.553

0.780

parenchyma-related

Kidney parenchymal (KT1)

Followup

0.968

0.389

0.518

Felzartamab

1.628

0.218

0.518

Followup:Felzartamab

0.605

0.552

0.552

Kidney parenchymal - no solute carriers (KT2)

Followup

0.859

0.432

0.518

Felzartamab

0.911

0.352

0.518

Followup:Felzartamab

1.038

0.365

0.518

1.6 SUMMARISE THE EFFECT OF TREATMENT FOR EACH SCORE ACROSS ALL TIME PERIODS


Here we present the effect of treatment (ΔΔ) for each score, along with the effect of time (Δ) for placebo and felzartmab arms.

flextable_pairwise

Table 1. Effect of felzartamab on molecular scores

Annotation

Score

Baseline - Week 24

Week 24 - Week 52

Baseline - Week 52

Δ Placebo
(N=10)

Δ Felzartamab
(N=10)

ΔΔ

ΔΔ FDR

Δ Placebo
(N=10)

Δ Felzartamab
(N=10)

ΔΔ

ΔΔ FDR

Δ Placebo
(N=10)

Δ Felzartamab
(N=10)

ΔΔ

ΔΔ FDR

cfDNA

Donor-derived cell-free DNA (dd-cfDNA, cp/mL)

0.50 (76.12)

-55.00 (30)

-55.50 (53.06)

2e-04

-3.50 (69.12)

14.50 (25.5)

18.00 (47.31)

0.027

-3.00 (51)

-40.50 (45.75)

-37.50 (48.38)

0.045

ABMR-related

ABMR classifier (ABMRProb)

0.10 (0.37)

-0.57 (0.33)

-0.67 (0.35)

1e-03

-0.26 (0.44)

0.54 (0.45)

0.79 (0.45)

1e-03

-0.15 (0.43)

-0.03 (0.42)

0.12 (0.43)

0.709

Glomerulitis classifier (g>0Prob)

-0.02 (0.23)

-0.31 (0.35)

-0.29 (0.29)

3e-03

-0.19 (0.34)

0.31 (0.41)

0.50 (0.37)

2e-04

-0.22 (0.3)

0.00 (0.38)

0.21 (0.34)

0.277

Peritubular capillaritis classifier (ptc>0Prob)

-0.12 (0.25)

-0.55 (0.27)

-0.43 (0.26)

3e-03

-0.10 (0.37)

0.37 (0.31)

0.47 (0.34)

3e-03

-0.22 (0.3)

-0.18 (0.35)

0.04 (0.33)

0.977

DSA-selective transcripts (DSAST)

-0.09 (0.22)

-0.37 (0.24)

-0.28 (0.23)

0.026

-0.16 (0.32)

0.16 (0.33)

0.32 (0.33)

0.026

-0.24 (0.33)

-0.20 (0.25)

0.04 (0.29)

0.988

NK cell burden (NKB)

0.02 (0.22)

-0.78 (0.49)

-0.80 (0.36)

0.014

-0.19 (0.4)

0.68 (0.59)

0.87 (0.49)

0.014

-0.17 (0.46)

-0.10 (0.62)

0.07 (0.54)

0.972

TCMR-related

TCMR classifier (TCMRProb)

-0.01 (0.02)

0.00 (0)

0.01 (0.01)

0.358

0.00 (0.02)

0.00 (0)

0.00 (0.01)

0.856

-0.01 (0.02)

0.00 (0)

0.01 (0.01)

0.358

Tubulitis classifier (t>1Prob)

-0.01 (0.06)

-0.02 (0.04)

-0.01 (0.05)

1.000

0.00 (0.04)

0.01 (0.02)

0.01 (0.03)

1.000

-0.01 (0.05)

-0.01 (0.03)

-0.01 (0.04)

1.000

Interstitial infiltrate classifier (i>1Prob)

-0.01 (0.04)

-0.02 (0.03)

-0.01 (0.04)

0.964

0.01 (0.04)

0.01 (0.04)

0.00 (0.04)

0.964

0.00 (0.05)

-0.01 (0.06)

-0.01 (0.05)

0.964

Cytotoxic T cell-associated (QCAT)

-0.20 (0.41)

-0.70 (0.46)

-0.50 (0.43)

0.425

0.10 (0.33)

0.54 (0.57)

0.44 (0.45)

0.425

-0.10 (0.46)

-0.16 (0.41)

-0.06 (0.44)

0.871

T cell burden (TCB)

-0.44 (0.71)

-0.79 (0.56)

-0.36 (0.64)

0.825

0.28 (0.57)

0.34 (0.6)

0.06 (0.58)

0.905

-0.16 (0.88)

-0.45 (0.43)

-0.29 (0.65)

0.825

macrophage-related

Alternatively activated macrophage (AMAT1)

0.08 (0.2)

-0.27 (0.45)

-0.35 (0.33)

0.128

0.09 (0.23)

0.11 (0.47)

0.02 (0.35)

0.812

0.18 (0.19)

-0.16 (0.32)

-0.33 (0.25)

0.128

Constitutive macrophage (QCMAT)

0.00 (0.15)

-0.24 (0.27)

-0.24 (0.21)

0.340

0.01 (0.12)

0.09 (0.27)

0.08 (0.2)

0.685

0.01 (0.11)

-0.15 (0.24)

-0.16 (0.18)

0.348

injury-related

Injury-repair associated (IRRAT30)

0.18 (0.52)

0.05 (0.4)

-0.13 (0.46)

0.353

0.05 (0.66)

-0.11 (0.72)

-0.16 (0.69)

0.353

0.23 (0.6)

-0.06 (0.48)

-0.29 (0.54)

0.175

Injury-repair induced, day 3 (IRITD3)

0.03 (0.19)

0.04 (0.12)

0.01 (0.15)

0.945

0.03 (0.18)

-0.18 (0.18)

-0.20 (0.18)

0.060

0.06 (0.19)

-0.13 (0.14)

-0.19 (0.16)

0.060

Injury-repair induced, day 5 (IRITD5)

-0.01 (0.17)

0.04 (0.14)

0.05 (0.15)

0.932

0.07 (0.2)

-0.17 (0.21)

-0.24 (0.21)

0.093

0.06 (0.18)

-0.13 (0.17)

-0.18 (0.18)

0.093

Fibrosis classifier (ci>1Prob)

-0.02 (0.29)

0.13 (0.32)

0.15 (0.3)

0.434

0.06 (0.25)

-0.08 (0.53)

-0.14 (0.39)

0.844

0.04 (0.25)

0.05 (0.37)

0.01 (0.31)

0.434

Atrophy classifier (ct>1Prob)

0.01 (0.33)

0.12 (0.31)

0.12 (0.32)

0.587

0.09 (0.27)

-0.06 (0.44)

-0.15 (0.36)

0.587

0.10 (0.32)

0.07 (0.41)

-0.03 (0.36)

0.587

parenchyma-related

Kidney parenchymal (KT1)

-0.05 (0.22)

-0.15 (0.23)

-0.09 (0.22)

0.731

-0.01 (0.11)

0.10 (0.22)

0.11 (0.17)

0.705

-0.06 (0.25)

-0.04 (0.18)

0.02 (0.22)

0.705

Kidney parenchymal - no solute carriers (KT2)

-0.09 (0.39)

-0.16 (0.37)

-0.07 (0.38)

0.577

-0.09 (0.2)

0.14 (0.38)

0.23 (0.29)

0.484

-0.18 (0.42)

-0.02 (0.31)

0.16 (0.37)

0.577

Grey shading denotes ANOVA interactive effect FDR < 0.05
FDR correction was carried out within each annotation grouping

1.7 PLOT THE RESULTS

# EXTRACT LEGEND FOR PLOTS ####
panel_legend <- plots_artANOVA %>%
    dplyr::filter(variable == "AMAT1") %>%
    pull(plot_violin) %>%
    ggpubr::get_legend() %>%
    ggpubr::as_ggplot() +
    theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))


# MOLECULAR ABMR PANELS ####
panel_violin_abmr <- plots_artANOVA %>%
    dplyr::filter(category %in% c("ABMR")) %>%
    pull(plot_violin) %>%
    wrap_plots(nrow = 1, ncol = 5) &
    theme(
        axis.text = element_text(size = 10, colour = "black"),
        legend.position = "none",
        plot.background = element_rect(fill = "grey95", colour = " white")
    )

panel_violin_abmr <- panel_violin_abmr %>%
    ggarrange(
        labels = "B",
        font.label = list(size = 25, face = "bold"),
        legend = "none"
    ) %>%
    ggpubr::annotate_figure(
        top = text_grob("Effect of felzartamab on molecular ABMR activity scores",
            face = "bold.italic",
            size = 25,
            hjust = 1.27
        )
    )


# MOLECULAR TCMR PANELS ####
panel_violin_tcmr <- plots_artANOVA %>%
    dplyr::filter(category %in% c("TCMR")) %>%
    pull(plot_violin) %>%
    wrap_plots(nrow = 1, ncol = 5) &
    theme(
        axis.text = element_text(size = 10, colour = "black"),
        legend.position = "none",
    )

panel_violin_tcmr <- panel_violin_tcmr %>%
    ggarrange(
        labels = "C",
        font.label = list(size = 25, face = "bold"),
        legend = "none"
    ) %>%
    ggpubr::annotate_figure(
        top = text_grob("Effect of felzartamab on molecular TCMR activity scores",
            face = "bold.italic",
            size = 25,
            hjust = 1.27
        )
    )


# COMBINED VIOLIN PANELS ####
panels_violin <- ggarrange(
    panel_violin_abmr,
    panel_violin_tcmr,
    nrow = 2,
    heights = c(1, 1)
)

ggarrange(
    panel_legend,
    panels_violin,
    nrow = 2,
    heights = c(0.125, 1)
) 

2 Analysing the effect of treatment on gene expression

2.1 INTRODUCTION

This analysis assesses how felzartamab therapy affects gene expression from biopsies from kidney allografts taken from patients with antibody mediated rejection. We use the Linear Models for Microarray Data (limma) pipeline for applying empirical Bayes mediated t-tests on individual genes. Specific contrasts to assessed the interaction between treatment and time are specified to estimate the effect of treatemnt on each gene.

References:

• Ritchie, M.E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).

2.2 DATA WRANGLING


First we need to load the necessary R packages and data.

# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse) 
library(flextable) 
library(officer)
# Bioconductor libraries
library(Biobase) 
library(limma) 
library(genefilter) 
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/data/affymap219.RData")

# Seed for reproducibility
set.seed(42)

# DEFINE THE SET ####
set00 <- set[, set$Patient %nin% c(15, 18)]


• We filter the gene expression data by to remove genes with minimal variance across the population (IQR filtering).

# IQR FILTER THE DATA ####
f1 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1)
if (!exists("selected")) {
    selected <- genefilter(set00, ff)
}
set01 <- set00[selected, ]

• We then retain a single probeset for each gene based on whichever has the highest mean expression in the population.

# KEEP UNIQUE GENES (keep probe with highest mean expression) ####
mean_exprs_by_probe <- set01 %>%
    exprs() %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
    mutate(mean_exprs = set01 %>%
        exprs() %>% rowMeans(), .after = Symb)

genes <- mean_exprs_by_probe %>%
    group_by(Symb) %>%
    dplyr::slice_max(mean_exprs) %>%
    dplyr::filter(Symb != "") %>%
    distinct(Symb, .keep_all = TRUE) %>%
    pull(AffyID)

set <- set01[featureNames(set01) %in% genes, ]

2.3 REMOVE GENES THAT DIFFER AT BASELINE BETWEEN PLACEBO AND FELZARTAMAB ARMS


The felzartamab trial was a phase II randomize control trial (RCT). While RCTs are the ‘gold standard’ for establishing causality (i.e., the effect of treatment), they are not perfect, particularly when the number of study participants is low. Thus, we pre-screened gene expression data to remove any genes that differed between placebo and felzartamab arm as baseline. This was done to reduce the potential for interpreting the effect of treatment on genes with dissimilar baseline profiles.

# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set$Felzartamab_Followup %>% droplevels()


# BLOCK DESIGN week24 - baseline ####
design_block01 <- model.matrix(~ 0 + Felzartamab_Followup)
contrast_block_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupBaseline_Felzartamab) -(Felzartamab_FollowupBaseline_Placebo)",
    levels = design_block01
)

# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_block_1 <- limma::lmFit(set, design_block01)
cfit_block_1 <- limma::contrasts.fit(fit_block_1, contrast_block_01)
ebayes_block_1 <- limma::eBayes(cfit_block_1)
tab_block_1 <- limma::topTable(ebayes_block_1, adjust = "fdr", sort.by = "p", number = "all")

# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means_baseline <- fit_block_1 %>%
    avearrays() %>%
    data.frame() %>%
    rownames_to_column("AffyID") %>%
    tibble() %>%
    mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
    rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
    dplyr::select(AffyID, contains("Baseline"))


The top 20 genes that differ at baseline (sorted by p-value) between placebo and felzartamab arms are summarized below.

(Note that all FDR > 0.05. This is a consequence of differential expression among with n = 10 by treatment group across a large selection of genes (5,267 in this case). Nonetheless, genes with p < 0.05 will be excluded from subsequent analyses).

# FORMAT TOPTABLES ####
baseline_deg <- tab_block_1 %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID") %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means_baseline, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb,
        Baseline_Placebo, Baseline_Felzartamab,
        t, logFC, P.Value, adj.P.Val,
    ) %>%
    dplyr::rename(p = P.Value, FDR = adj.P.Val)

baseline_deg %>%
    dplyr::slice_min(p, n = 20) %>%
    mutate_at(vars(t, logFC, FDR), ~ formatC(.x, format = "f", digits = 2)) %>%
    mutate_at(vars(p), ~ formatC(.x, format = "f", digits = 5)) %>%
    flextable::flextable()

AffyID

Symb

Baseline_Placebo

Baseline_Felzartamab

t

logFC

p

FDR

11727331_at

SPAG4

90

50

-4.23

-0.85

0.00008

0.33

11763725_a_at

HNRNPU

366

665

4.10

0.86

0.00013

0.33

11718286_a_at

UBL3

74

142

3.66

0.94

0.00053

0.45

11732590_x_at

TRIM5

42

63

3.65

0.59

0.00055

0.45

11725671_a_at

IPO7

42

66

3.54

0.64

0.00079

0.45

11741431_a_at

KCNJ16

203

311

3.49

0.61

0.00091

0.45

11726915_at

SLCO4C1

87

176

3.48

1.02

0.00096

0.45

11723109_a_at

GPBP1

153

263

3.46

0.78

0.00099

0.45

11739744_x_at

MXI1

123

191

3.45

0.64

0.00102

0.45

11720890_a_at

AGL

39

63

3.40

0.70

0.00121

0.45

11746431_a_at

FBXO11

28

45

3.38

0.68

0.00127

0.45

11751841_a_at

TNFRSF17

56

20

-3.36

-1.51

0.00134

0.45

11744367_a_at

ZMYND11

30

49

3.33

0.68

0.00150

0.45

11748857_a_at

AK3

371

578

3.32

0.64

0.00155

0.45

11727307_x_at

ZNF320

46

66

3.31

0.51

0.00160

0.45

11747648_a_at

TCAIM

96

152

3.29

0.66

0.00168

0.45

11728966_s_at

ZNF28

29

67

3.27

1.21

0.00181

0.45

11723387_a_at

PTAR1

149

261

3.25

0.81

0.00187

0.45

11739771_a_at

OXR1

31

51

3.24

0.73

0.00194

0.45

11730038_at

GNG13

44

28

-3.24

-0.63

0.00196

0.45

2.4 ESTIMATING TREATMENT EFFECT ON GENE EXPRESSION


We now carry out differential expression analysis to estimate the treatment effect on individual genes.

# DEFINE GENES SIMILAR AT BASELINE ####
genes_baseline <- baseline_deg %>%
    dplyr::filter(p > 0.05) %>%
    pull(AffyID)

# WRANGLE THE MEAN EXPRESSION DATA ####
mean_exprs_by_probe <- set01 %>%
    exprs() %>%
    as_tibble(rownames = "AffyID") %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
    mutate(
        mean_exprs = set01 %>%
            exprs() %>%
            rowMeans(),
        .after = Symb
    )

genes <- mean_exprs_by_probe %>%
    group_by(Symb) %>%
    dplyr::slice_max(mean_exprs) %>%
    dplyr::filter(Symb != "", AffyID %in% genes_baseline) %>%
    distinct(Symb, .keep_all = TRUE) %>%
    pull(AffyID)

set02 <- set01[featureNames(set01) %in% genes, ]


We need to define the model structure for the contrasts we intend to make. Here we will define 3 major contrast groups; 1) effect of time in placebo group, 2) effect of time in felzartamab group, and 3) the interactive effect of time and treatment. We do this for each time window (baseline to week 24, week 24 to week 52, and baseline to week 52)

# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set02$Felzartamab_Followup %>% droplevels()
cortex <- set02$Cortexprob


# DESIGN ####
design <- model.matrix(~ 0 + Felzartamab_Followup + cortex)


# CONTRAST DESIGN week24 - baseline ####
contrast_interaction_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)
contrast_felzartamab_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
    levels = design
)
contrast_placebo_01 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)


# CONTRAST DESIGN week52 - week24 ####
contrast_interaction_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
    levels = design
)
contrast_felzartamab_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2",
    levels = design
)
contrast_placebo_02 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
    levels = design
)


# CONTRAST DESIGN week52 - baseline####
contrast_interaction_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)
contrast_felzartamab_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
    levels = design
)
contrast_placebo_03 <- makeContrasts(
    "x =  (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
    levels = design
)


We now fit the models we’ve specified

# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_interaction_1 <- limma::lmFit(set02, design)
cfit_interaction_1 <- limma::contrasts.fit(fit_interaction_1, contrast_interaction_01)
ebayes_interaction_1 <- limma::eBayes(cfit_interaction_1)
tab_interaction_1 <- limma::topTable(ebayes_interaction_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_1$stdev.unscaled * sqrt(ebayes_interaction_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_1 <- limma::lmFit(set02, design)
cfit_felzartamab_1 <- limma::contrasts.fit(fit_felzartamab_1, contrast_felzartamab_01)
ebayes_felzartamab_1 <- limma::eBayes(cfit_felzartamab_1)
tab_felzartamab_1 <- limma::topTable(ebayes_felzartamab_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_1$stdev.unscaled * sqrt(ebayes_felzartamab_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_1 <- limma::lmFit(set02, design)
cfit_placebo_1 <- limma::contrasts.fit(fit_placebo_1, contrast_placebo_01)
ebayes_placebo_1 <- limma::eBayes(cfit_placebo_1)
tab_placebo_1 <- limma::topTable(ebayes_placebo_1, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_1$stdev.unscaled * sqrt(ebayes_placebo_1$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")


# FIT BLOCK week52 - week24 LIMMA MODEL ####
fit_interaction_2 <- limma::lmFit(set02, design)
cfit_interaction_2 <- limma::contrasts.fit(fit_interaction_2, contrast_interaction_02)
ebayes_interaction_2 <- limma::eBayes(cfit_interaction_2)
tab_interaction_2 <- limma::topTable(ebayes_interaction_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_2$stdev.unscaled * sqrt(ebayes_interaction_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_2 <- limma::lmFit(set02, design)
cfit_felzartamab_2 <- limma::contrasts.fit(fit_felzartamab_2, contrast_felzartamab_02)
ebayes_felzartamab_2 <- limma::eBayes(cfit_felzartamab_2)
tab_felzartamab_2 <- limma::topTable(ebayes_felzartamab_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_2$stdev.unscaled * sqrt(ebayes_felzartamab_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_2 <- limma::lmFit(set02, design)
cfit_placebo_2 <- limma::contrasts.fit(fit_placebo_2, contrast_placebo_02)
ebayes_placebo_2 <- limma::eBayes(cfit_placebo_2)
tab_placebo_2 <- limma::topTable(ebayes_placebo_2, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_2$stdev.unscaled * sqrt(ebayes_placebo_2$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")


# FIT BLOCK week52 - baseline LIMMA MODEL ####
fit_interaction_3 <- limma::lmFit(set02, design)
cfit_interaction_3 <- limma::contrasts.fit(fit_interaction_3, contrast_interaction_03)
ebayes_interaction_3 <- limma::eBayes(cfit_interaction_3)
tab_interaction_3 <- limma::topTable(ebayes_interaction_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_interaction_3$stdev.unscaled * sqrt(ebayes_interaction_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_felzartamab_3 <- limma::lmFit(set02, design)
cfit_felzartamab_3 <- limma::contrasts.fit(fit_felzartamab_3, contrast_felzartamab_03)
ebayes_felzartamab_3 <- limma::eBayes(cfit_felzartamab_3)
tab_felzartamab_3 <- limma::topTable(ebayes_felzartamab_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_felzartamab_3$stdev.unscaled * sqrt(ebayes_felzartamab_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

fit_placebo_3 <- limma::lmFit(set02, design)
cfit_placebo_3 <- limma::contrasts.fit(fit_placebo_3, contrast_placebo_03)
ebayes_placebo_3 <- limma::eBayes(cfit_placebo_3)
tab_placebo_3 <- limma::topTable(ebayes_placebo_3, adjust = "fdr", sort.by = "p", number = "all") %>%
    as_tibble(rownames = "AffyID") %>%
    mutate(se = (ebayes_placebo_3$stdev.unscaled * sqrt(ebayes_placebo_3$s2.post)) %>% as.vector(), .after = logFC) %>%
    right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")

2.5 ORAGNIZE THE RESULTS


The results are now organized for tabulation

# MERGE BLOCKS ####
tab_block_1 <- tab_interaction_1 %>%
    left_join(
        tab_felzartamab_1 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_1 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)

tab_block_2 <- tab_interaction_2 %>%
    left_join(
        tab_felzartamab_2 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_2 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)

tab_block_3 <- tab_interaction_3 %>%
    left_join(
        tab_felzartamab_3 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(flogFC = logFC),
        by = "AffyID"
    ) %>%
    left_join(
        tab_placebo_3 %>%
            dplyr::select(AffyID, logFC) %>%
            rename(plogFC = logFC),
        by = "AffyID"
    ) %>%
    relocate(plogFC, flogFC, .before = logFC) %>%
    arrange(P.Value)


# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means <- fit_interaction_1 %>%
    avearrays() %>%
    as_tibble(rownames = "AffyID") %>%
    mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
    rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
    dplyr::select(-contains("Week12"), -any_of(c("cortex")))


# FORMAT TOPTABLES ####
table_interaction_1 <- tab_block_1 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

table_interaction_2 <- tab_block_2 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

table_interaction_3 <- tab_block_3 %>%
    arrange(P.Value) %>%
    mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
    tibble() %>%
    left_join(means, by = "AffyID") %>%
    dplyr::select(
        AffyID, Symb, Gene, PBT,
        t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
        all_of(colnames(means))
    ) %>%
    dplyr::rename(
        p = P.Value,
        FDR = adj.P.Val
    )

deg <- tibble(
    design = c(
        "Baseline_vs_Week24",
        "Week24_vs_Week52",
        "Baseline_vs_Week52"
    ),
    table = list(
        table_interaction_1,
        table_interaction_2,
        table_interaction_3
    )
)

2.6 SUMMARIZE THE RESULTS

# FORMAT TABLES TO MAKE FLEXTABLES ####
table_deg <- deg %>%
    expand_grid(direction = c("all", "increased", "decreased")) %>%
    relocate(direction, .after = design) %>%
    mutate(
        gene_tables = pmap(
            list(direction, table),
            function(direction, table) {
                if (direction == "increased") {
                    table <- table %>%
                        dplyr::filter(logFC > 0) %>%
                        arrange(p)
                } else if (direction == "decreased") {
                    table <- table %>%
                        dplyr::filter(logFC < 0) %>%
                        arrange(p)
                }
                table %>%
                    dplyr::select(-t, -se, -FDR, -contains("AffyID"), -contains("MMDx")) %>%
                    dplyr::slice(1:20) %>%
                    mutate(
                        Gene = Gene %>% str_remove("///.*"),
                        plogFC = plogFC %>% round(2),
                        flogFC = flogFC %>% round(2),
                        logFC = logFC %>% round(2),
                        p = case_when(
                            p < 0.0001 ~ p %>% formatC(digits = 0, format = "e"),
                            TRUE ~ p %>% formatC(digits = 4, format = "f")
                        )
                    )
            }
        )
    )


# GLOBAL PARAMETERS FOR FLEXTABLES ####
header1 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    rep("Differential expression", 4),
    rep("Mean expression by group", 6)
)
header2 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    "\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
    "\u394\u394\nlogFC", "\u394\u394\nP",
    rep("Placebo", 3), rep("Felzartamab", 3)
)
header3 <- c(
    "Gene\nsymbol", "Gene", "PBT",
    "\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
    "\u394\u394\nlogFC", "\u394\u394\nP",
    "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)", "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)"
)

# MAKE FORMATTED FLEXTABLES ####
flextables <- table_deg %>%
    mutate(
        flextable = pmap(
            list(design, gene_tables),
            function(design, gene_tables) {
                if (design == "Baseline_vs_Week24") {
                    title <- paste("Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                } else if (design == "Week24_vs_Week52") {
                    title <- paste("Table i. Top 20 differentially expressed genes between week24 and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                } else if (design == "Baseline_vs_Week52") {
                    title <- paste("Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
                }
                gene_tables %>%
                    arrange(logFC) %>%
                    flextable::flextable() %>%
                    flextable::delete_part("header") %>%
                    flextable::add_header_row(top = TRUE, values = header3) %>%
                    flextable::add_header_row(top = TRUE, values = header2) %>%
                    flextable::add_header_row(top = TRUE, values = header1) %>%
                    flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
                    flextable::merge_v(part = "header") %>%
                    flextable::merge_h(part = "header") %>%
                    flextable::border_remove() %>%
                    flextable::border(part = "header", border = officer::fp_border()) %>%
                    flextable::border(part = "body", border = officer::fp_border()) %>%
                    flextable::align(align = "center") %>%
                    flextable::align(align = "center", part = "header") %>%
                    flextable::font(fontname = "Arial", part = "all") %>%
                    flextable::fontsize(size = 7, part = "all") %>%
                    flextable::fontsize(size = 7, part = "footer") %>%
                    flextable::fontsize(i = 1, size = 12, part = "header") %>%
                    flextable::bold(part = "header") %>%
                    flextable::bg(bg = "white", part = "all") %>%
                    flextable::padding(padding = 0, part = "all") %>%
                    flextable::autofit()
            }
        )
    )


Genes suppressed by felzartamab in week 24 compared to baseline
(Note that all ΔΔFDR > 0.05. This is a consequence of differential expression among with n = 10 by treatment group across a large selection of genes (5,267 in this case). Thus, treatment effects on individual genes should not be interpreted. Instead, geneset enrichment analysis (GSEA) will be used to assess bulk trends).

flextables %>%
    dplyr::filter(design == "Baseline_vs_Week24", direction == "decreased") %>%
    pull(flextable) %>%
    pluck(1)

Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)

Gene
symbol

Gene

PBT

Differential expression

Mean expression by group

Δ
placebo
logFC

Δ
felzartamab
logFC

ΔΔ
logFC

ΔΔ
P

Placebo

Felzartamab

Baseline
(N=10)

Week24
(N=10)

Week52
(N=10)

Baseline
(N=10)

Week24
(N=10)

Week52
(N=10)

CXCL11

chemokine (C-X-C motif) ligand 11

ABMR-RAT,GRIT3,RAT,Rej-RAT

-0.16

-1.01

-0.85

0.0176

494

395

372

617

153

317

CXCL9

chemokine (C-X-C motif) ligand 9

ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT

-0.14

-0.94

-0.79

0.0290

1,344

1,102

1,384

1,684

459

1,036

CXCL10

chemokine (C-X-C motif) ligand 10

ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT

-0.21

-0.98

-0.77

0.0228

844

632

720

1,063

274

615

IDO1

indoleamine 2,3-dioxygenase 1

ABMR-RAT,GRIT3,RAT,Rej-RAT

-0.09

-0.77

-0.68

0.0103

338

298

296

409

140

251

LYZ

lysozyme

IRITD5,QCMAT

-0.01

-0.67

-0.66

0.0386

1,676

1,658

1,884

2,244

884

1,184

FCGR3A

Fc fragment of IgG, low affinity IIIa, receptor (CD16a)

ABMR-RAT,GRIT2,IRRAT950,RAT,Rej-RAT

0.06

-0.58

-0.64

0.0167

550

602

613

571

256

352

GZMB

granzyme B

ABMR-RAT,QCAT,RAT,Rej-RAT

-0.01

-0.62

-0.61

0.0019

149

147

130

143

61

114

GBP1

guanylate binding protein 1, interferon-inducible

ABMR-RAT,GRIT3,RAT,Rej-RAT

-0.03

-0.61

-0.59

0.0096

678

655

643

909

389

556

GNLY

granulysin

ABMR-RAT,DSAST,QCAT,RAT,Rej-RAT

-0.18

-0.75

-0.57

0.0285

192

149

134

165

58

100

SH2D1B

SH2 domain containing 1B

ABMR-RAT,DSAST,NKB,RAT

0.07

-0.44

-0.51

0.0010

38

42

34

40

22

33

FCGR1A

Fc fragment of IgG, high affinity Ia, receptor (CD64)

GRIT2,IRRAT950,RAT,TCMR-RAT

0.09

-0.41

-0.50

0.0304

131

148

161

152

86

86

KLRD1

killer cell lectin-like receptor subfamily D, member 1

ABMR-RAT,RAT,Rej-RAT

-0.06

-0.49

-0.43

0.0295

110

101

96

126

63

102

EPB41L3

erythrocyte membrane protein band 4.1-like 3

0.22

-0.20

-0.42

0.0388

209

286

258

339

258

283

LOC339059

uncharacterized LOC339059

0.26

-0.15

-0.41

0.0203

51

73

55

66

53

59

APOL1

apolipoprotein L1

ABMR-RAT,GRIT3,RAT,Rej-RAT

0.07

-0.32

-0.39

0.0205

626

687

744

754

483

636

LYPD5

LY6/PLAUR domain containing 5

ABMR-RAT,RAT,Rej-RAT

0.04

-0.33

-0.37

0.0255

27

28

27

30

19

26

TRAV12-2

T cell receptor alpha variable 12-2

0.25

-0.10

-0.35

0.0173

39

55

39

49

43

43

IL18RAP

interleukin 18 receptor accessory protein

0.02

-0.32

-0.34

0.0233

38

39

41

42

27

33

PRF1

perforin 1 (pore forming protein)

ABMR-RAT,QCAT,RAT,Rej-RAT

-0.02

-0.35

-0.33

0.0382

152

148

137

145

89

123

PMM2

phosphomannomutase 2

0.20

-0.07

-0.27

0.0380

38

50

41

46

42

49


Genes suppressed by felzartamab in week 52 compared to baseline
(Note that all ΔΔFDR > 0.05. This is a consequence of differential expression among with n = 10 by treatment group across a large selection of genes (5,267 in this case). Thus, treatment effects on individual genes should not be interpreted. Instead, geneset enrichment analysis (GSEA) will be used to assess bulk trends).

flextables %>%
    dplyr::filter(design == "Baseline_vs_Week52", direction == "decreased") %>%
    pull(flextable) %>%
    pluck(1)

Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)

Gene
symbol

Gene

PBT

Differential expression

Mean expression by group

Δ
placebo
logFC

Δ
felzartamab
logFC

ΔΔ
logFC

ΔΔ
P

Placebo

Felzartamab

Baseline
(N=10)

Week24
(N=10)

Week52
(N=10)

Baseline
(N=10)

Week24
(N=10)

Week52
(N=10)

HNRNPK

heterogeneous nuclear ribonucleoprotein K

0.28

-0.38

-0.66

0.0055

818

1,044

1,203

1,096

1,059

649

COL1A1

collagen, type I, alpha 1

COLA1356,FICOL,IRITD5,LivGST_UP

0.38

-0.26

-0.64

0.0075

310

333

526

329

320

230

NFIX

nuclear factor I/X (CCAAT-binding transcription factor)

CT1

0.42

-0.15

-0.57

0.0012

162

179

289

221

167

180

RHOA

ras homolog family member A

IRRAT950

0.28

-0.29

-0.57

0.0049

691

803

1,017

923

809

620

COL1A2

collagen, type I, alpha 2

FICOL,IRITD5,LivGST_UP

0.26

-0.31

-0.57

0.0118

1,238

1,152

1,769

1,534

1,194

997

ITGB3

integrin beta 3

IRRAT30,IRRAT950

0.34

-0.21

-0.55

0.0100

105

137

169

110

103

82

LOC100129518

uncharacterized LOC100129518

GRIT3,IRRAT950

0.25

-0.22

-0.47

0.0050

429

546

605

527

461

390

ACTB

actin, beta

IRRAT950

0.15

-0.31

-0.47

0.0077

6,451

7,352

7,993

7,274

5,662

4,720

SPARC

secreted protein, acidic, cysteine-rich (osteonectin)

IRITD3

0.16

-0.30

-0.46

0.0100

2,411

2,532

3,011

2,343

2,094

1,550

DNAJA4

DnaJ (Hsp40) homolog, subfamily A, member 4

IRRAT950

0.27

-0.14

-0.42

0.0095

49

68

72

68

61

56

SLFN5

schlafen family member 5

0.19

-0.20

-0.40

0.0050

39

44

51

49

38

37

CYR61

cysteine-rich, angiogenic inducer, 61

IRITD3,IRRAT950

0.14

-0.26

-0.40

0.0074

271

268

330

338

234

236

SDCBP

syndecan binding protein

0.24

-0.14

-0.38

0.0120

1,135

1,404

1,582

1,478

1,358

1,224

SBSPON

somatomedin B and thrombospondin type 1 domain containing

0.31

-0.08

-0.38

0.0123

217

250

332

188

249

169

PAFAH1B2

platelet-activating factor acetylhydrolase 1b, catalytic subunit 2 (30kDa)

0.23

-0.13

-0.36

0.0099

309

372

423

386

396

323

TMEM2

transmembrane protein 2

0.10

-0.26

-0.36

0.0118

270

307

311

318

259

222

CALM2

calmodulin 2 (phosphorylase kinase, delta)

CT1,IRITD3

0.12

-0.25

-0.36

0.0138

1,304

1,427

1,538

1,623

1,342

1,155

WNK1

WNK lysine deficient protein kinase 1

0.26

-0.08

-0.34

0.0104

703

884

1,012

894

996

799

CEBPB

CCAAT/enhancer binding protein (C/EBP), beta

IRITD3

0.11

-0.22

-0.33

0.0123

419

452

488

417

333

307

STAT6

signal transducer and activator of transcription 6, interleukin-4 induced

0.09

-0.20

-0.30

0.0144

671

710

765

746

616

563